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We study 2D fronts propagating up a co-moving reaction rate gradient in finite number reaction- 
diffusion systems. We sliow that in a 2D rectangular channel, planar solutions to the deterministic 
mean-field equation are stable with respect to deviations from planarity. We argue that planar 
fronts in the corresponding stochastic system, on the other hand, are unstable if the channel width 
exceeds a critical value. Furthermore, the velocity of the stochastic fronts is shown to depend on 
the channel width in a simple and interesting way, in contrast to fronts in the deterministic MFE. 
Thus, fiuctuations alter the behavior of these fronts in an essential way. These affects are shown to 
be partially captured by introducing a density cutoff in the reaction rate. Some of the predictions 
of the cutoff mean-field approach are shown to be in quantitative accord with the stochastic results. 
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I. INTRODUCTION 

Several well known processes in spatially extended systems exhibit fronts that propagate through space. Most 
of these processes that have been considered to date occur in media in which the governing dynamics are spatially 
uniform. Recently, however, some interesting findings have been made concerning fronts propagating in systems with 
spatially heterogeneous dynamics. In particular, the simple infection model A -f i? ^ 2A on a lattice with equal 
hopping rates and a linear reaction rate gradient has been studied. Two versions of this system have been examined 
in some detail PjQ ■ oiic in which the gradient is defined with respect to the medium itself (the "absolute gradient"), 
and another in which the gradient is defined relative the the front's interface and travels along with the front (the 
"quasi-static gradient"). One can imagine numerous systems that can be described by the absolute gradient, e.g. a 
chemical reaction occurring in a temperature gradient. The quasi-static gradient is more analytically tractable and 
also arises naturally in models of biological evolution0,U. 

The usual way to analytically study a system with a propagating front, such as the infection model mentioned 
above, is within a mean field, reaction-diffusion framework. The simplest MP analog to our infection model is the 
usual Fisher equation l5|, with a spatially varying reaction rate: 

d(j)/dt = DV'^(j) + r{x)(l){l-(l)) (1) 

For our simple infection model, EqQ] (the " naive MFE" ) fails to to capture many of the qualitative aspects of the 
stochastic problem with either absolute of quasi-static gradients. These failures, as well as many other many other 
issues involving the MF description of similar front propagation problems, are largely remedied by introducing a cutoff 
factor in the reaction term: 

5(/)/5f = DV2<^-l-r(x)0(l-0)6'(0-0c) (2) 
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This term causes the reaction rate to abruptly drop to zero in regions far into the front's leading edge where (j) 
drops below a critical level <f>c, and is meant to mimic the leading order effect of finite number fluctuations in the 
stochastic process. In other words, the discrete nature of individual particles implies that a sufficiently small density 
field (j) < 4>c ^ 1/N corresponds, in an average sense, to zero particles present and thus zero reaction rate. 

In this study, we turn our attention to the quasi-static gradient in two dimensions and ask how finite number 
fluctuations and the related cutoff approach affect the stability of planar fronts propagating in a rectangular channel. 
In what follows we will see that in contrast to the predictions of the naive MFE, the results of stochastic simulations 
point to unstable planar fronts. Furthermore, once again the cutoff term will rescue the effectiveness of the mean 
field description of the stochastic process. We first study the cutoff mean-field equations, both numerically and 
analytically, showing the instability. We then turn to the stochastic model, demonstrating the instability there as 
well. An appendix contains details about the numerics. 

II. MEAN FIELD STABILITY CALCULATION 

The full equation of motion governing the quasi-static gradient in the MF cutoff framework is Eq. [21 with 

= Na/N 

r{x) = max(rmi„, To + a{x ~ x)) 

x = lJ (j)(x, y)dxdy 

0c = k/N 

Here x is the direction parallel to the channel's long axis, N is the equilibrium number of A particles per lattice site, 
and k is some 0(1) fitting parameter, x serves to define the interface position of the front by essentially comparing 
the front's profile to a step function, and h is the cross-channel width. Tmin merely serves to keep the front stable and 
plays no important role, as its effects are only felt far from the leading edge. In the context of biological evolution, 
represents the fraction of individuals in a population with a given fitness x. If the size of the population is fixed, the 
growth rate of individuals with a particular fitness is proportional to x — x. The diffusion term represents "genetic 
drift" , and the dynamics of the system corresponds to the population evolving towards greater mean fitness. 

Direct numerical integration of the time dependent Eq. [21 shows that planar fronts are in fact unstable. For a 
sufficiently wide channel, perturbed planar fronts develop into long, though finite, fingers whose length increases with 
increasing channel width. An example of such a finger is shown in Fig. ^ We see that there is a deep narrow "notch" 
on the trailing side of the finger, so that the width of the interface is much greater here than for the rest of the finger. 
Defining the finger length by / [0(a;,O) — </)(a;, 6)] dx, the data for finger length versus channel width is presented in 
Fig. [21 We now turn toward an analytic understanding of this result. 

Due to the translational invariance of the system, it is natural to investigate first steady-state propagating planar 
front solutions. Plugging into Eq. ^ the traveling wave form, (j)Q{x,y,t) = 4>q{x — vt), with velocity v, we obtain 

Dc^'^ + vc^'^ + r(z)0o(l - Me U^ - ^) = (3) 

in terms of the comoving coordinate z ^ x — vt. A quick analysis of the linearized version of Eq. ||2Il provides insight 
into the role of the cutoff. As z — > — oo, (pQ —f 1. Linearizing around 4>o = 1; we find two exponential solutions, but 
one must be discarded since it decreases with increasing z. Similarly, as z ^ oo, (jjo ~^ 0. Linearizing around 0o — 
0, in the region past the cutoff, once again we find only one acceptable, decaying solution. This leaves our solution 
with a total of two undetermined constants. Fixing translational invariance reduces this number to one. Requiring 
continuity of 4>o at the cutoff determines the remaining coefficient, and continuity of c/jq determines the velocity. Thus, 
mathematically, the cutoff fixes the velocity by overdetermining the boundary conditions, i.e. converting Eq. Q into 
an eigenvalue problem. An analysis for large N yields the leading order result [l|: 

v^[2AD^a\ii{N/k)]^^^ (4) 

In the limit k/N ^ we regain the naive MF approach, in which u — > oo. Thus the naive MF and the cutoff MF 
predict qualitatively different results with respect to velocity. Not surprisingly, stochastic fronts in fact approach a 
(finite) steady-state velocity, in line with the cutoff MF. 

Turning now to 2D fronts, we wish to study the linear stability of the planar solution to transverse perturbations. 
We write 4>{x,y,t) — (t>o{z) + ^'i^iVj't) ^-i^d linearize Eq. Q with respect to 4>. The invariance of the system with 
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FIG. 1: Snapshots of the growing finger for the cutoff MFE compared to that for the stochastic model. The parameters are 
L» = 1, ro = 6, Q = 0.3. a)The cutoff MFE with k/N = 8.7 • 10"^ b) the stochastic model with N = 90801, c) the stochastic 
model with A^ = 2881. 



respect to translations in time and the transverse spatial direction y implies (f>{z,y,t) = e"*e''^*'77(z). The governing 
equation for r]{z) is then: 



with 



Dtj" + vrl + 7?r(z)[(l - 20„)0((^o - kjN) + 0o(l - '/'o)'5(0o " fc/A^)] = ^^ 



n = Dq^ 



(5) 



(6) 



The delta function arises from differentiating the step function and is due to the shift in Zcut caused by the perturbation 
(j>. We have assumed here that g 7^ so that J (j){z, y, t)dy — 0. The case q = has to be treated separately, but 
in any case the least stable mode should be the translation mode with 17 = 0. Notice that Eq. IjHJl implies a simple 
stabilizing quadratic dependence of the growth rate lo on q. Thus the least stable mode is that with the smallest 
non-zero q, which, assuming a zero-flux sidewall boundary condition, is qmin = t^ /h. This implies a minimum channel 
width b* below which even the longest wavelength mode has too much curvature for any instability to exist: 



6* =7r 



D 



^r. 



(J) 



where 0,max is the largest eigenvalue of the stability operator, Eq. jSJ. 

Like the steady-state problem, insight can be gained into Eq. |SJl by considering the boundary conditions at 
z -^ ±cx3. We require that 77 ^ as 2: — > — cx). As 0o ^ 1 in this region, we find two exponential solutions for ry: 
one growing with increasing z and the other decaying. The latter must of course be excluded. If we perform the 
same procedure past the cutoff, as z — > 00, we find two decaying modes. However, one of the modes decays more 
slowly than the steady-state solution and thus dominates it for sufficiently large z. This is unacceptable behavior for 
a perturbation and therefore this solution is discarded. Thus, our solution has two arbitrary constants, one of which 
may be chosen arbitrarily since Eq. (O is linear in r]. The remaining constant is fixed by requiring continuity of r] at 
the cutoff. Matching ry' at the cutoff determines the eigenvalues il. Thus, once again the cutoff has played a central 
role in determining the problem's interesting quantities. 
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FIG. 2: Finger length (as defined in the text) versus channel width for the cutoff MFE. The parameters are: D — 1, ro — 6, 
a = 0.3, k/N = 8.7-10-^. 

As with the steady-state problem, we can make analytic progress in the limit of large N. In this case, the cutoff is 
at large z, in the region where (po is small. If we consider Eq. (O in the region where (pQ ^ I and z < Zcut, and fix 
the translation invariance by setting z = for the unperturbed state, we obtain 



Drj" + vrj' + ry(ro + az) ~ Clrj 
Up to a similarity transformation, this is the Airy equation, with the general solution 
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(8) 



(9) 



with 



_ v^/4D-ro+n 
a 



'5 -(f) 

We argue that the Bi term must vanish by considering the large v limit of Eq. lO and matching onto Eq. . As 
shown in [2I, in the large v limit, the diffusion term in Eq. Q can be ignored, and the solution in the region where 
(/)o < 1 is 



Tj ^ e 



-i[(n)-n)z+iQz^] 



(10) 



Expansions of Ai and Bi for large argument show that the diffusionless result Eq. (|10ll matches onto Eq. © only 
if the Bi term is absent. The constant A may be arbitrarily set to unity since the problem is linear. Thus we have 
for z < Zcut 



■q = e 2D Ai (^-7^) 



We have to match this result to the simple exponential solution for z > z^ut- Thus, 



V — Zr 



e-"^Ai r ^""M =Cel ^° V*'v-' "^ yj (11) 

The derivative of r\ must also match properly at the cutoff. Looking back to Eq. © , we see that the delta function 
term causes a discontinuity in rl at Z(^ut- 

, _, r(ze.O#(l-fc/iV) _ r{z,ut), , 

'^"0'^' ^ie/t- ^ ^,(^^^^) - ^ (1 fc/A/j (12) 

Computing the derivatives in (12) and dividing by (11), we obtain 

— ^l + 4Dn/v^ - -(ro + az™,)(l " k/N) = -—^-^^ (13) 

This equation determines Q if the quantities v and Zcut are known. Now, for large v the LHS of (12) is also large. For 
the RHS to balance it, the Airy function in the denominator must be small. Thus, ~^°"* ~ ^ where ^o ~ —2.3381 is 
the first zero of the Airy function. For the value of the cutoff, we quote another result from [l| obtained by matching 
the linearized steady-state equation at the cutoff: 

Zcut « ""'^^^ ~ '^" - ^oS - 2D/V (14) 

a 

Plugging this expression into Eq. (|13() . expanding around ^Oj and dropping higher order terms, we obtain the leading 
order result valid for large v. 

n='-^ (15) 

V 

Here we see the interesting result that plane fronts become stable in the limit as w -^ oo, i.e. N -^ oo, i.e. the cutoff 
disappears. Eq. H15|l can be interpreted as saying the that Q. is proportional to the ratio of the diffusive length scale 
(D/v) to the length scale over which the rate changes appreciably (1/a). Thus, heuristically incorporating the effects 
of finite number ffuctuations qualitatively changed the system's stability properties by limiting the front's velocity, 
which in turn makes the diffusive length scale finite. This fluctuation induced instability is similar to that in 6^, 
where it was found that that a coupled reaction diffusion system with no reaction gradient, but with unequal diffusion 
coefficients, is unstable with a cutoff but stable without one. Furthermore, Eq. I|15l) shows that the fronts become 
stable as a — > 0, for any value of N . This is analogous to the stability shown in 6] in the case of equal diffusion 
coefficients. 

Thus, once again, the presence of the cutoff qualitatively changes the simple mean field predictions. If the cutoff 
approach is to mimic the effects of finite number particle fluctuations, we should expect to see some analog of this 
front instability in the stochastic, discrete infection model discussed earlier, to which we now return. 

III. STOCHASTIC SYSTEM 

We ran simulations in which the lower rectangular portion of the channel was initially populated with N type A 
particles per lattice site and the upper rectangular portion of the channel which was populated by N type B particles 
per lattice site. During each time step, a binomially distributed random number of particles hops to adjacent sites. 
Furthermore, A particles probabilistically cause some B particles to change into A particles . The reaction probability 
and hopping rates were chosen so that the discretized, stochastic equation for ATV^ reduces to EqQ] (with the quasi- 
static form for r{x)) when the expectation value is taken in the small time, small lattice spacing limit. In particular, 
for the hopping probability we took Phop ^ D-p-, where dt is the simulation timestep and I is the lattice spacing. The 
number of particles reacting during each timestep was chosen as a binomially distributed random variable with mean 
1 _ (1 _ lM^-^Nu drawn Na times. 

Results of simulations depended crucially on the width of the channel b and, to a lesser extent, the equilibrium 
number of particles per lattice site A'^. Initially uniform planar fronts were roughened by random fluctuations. These 
fluctuations, especially those occurring near the interface, often developed into a more organized, longer wavelength 
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FIG. 3: (Color online)The circles represent the exact numeric solution of Eq. 0. The solid line is the exact numeric solution 
of Eq. 1131 which is itself a large v approximation. Exact numerically generated values of A'^ and Zcut were used to generate 
this approximation. The dashed curve is the analytic approximation Eq. p5|l . The parameters are D — tq = 1, a — 0.3. 



deviations from planarity. Naturally, these fluctuations occurred more frequently in runs with smaller N. For b <b*, 
diffusion tended to flatten these deviations back to planarity over a characteristic time scale tdev Although tdev is 
difficult to measure quantitatively, it is qualitatively clear that it increased with increasing b. Not surprisingly, as 
tdev and b increased, so did the amplitude of the associated deviations. The divergence of tdev corresponds to the 
aforementioned phase transition predicted by the cutoff MF approach in Eq. [7| For b ^ b* tdev did in fact appear 
to diverge, resulting in patterns that resemble fingers, rather that random fluctuations. As b is increased further, the 
finger-like nature of the pattern becomes increasingly pronounced. This trend is of course clearly the larger N is. In 
Fig. n we present the stochastic finger for two different values of N, alongside the pattern produced by the cutoff 
MF. The larger N corresponds to the value of the cutoff chosen for the MF simulation. The overall similarity between 
these two figures is clear. What is also striking is the relatively large effect of noise, even though N is quite large, 
pointing out once more the crucial importance of fluctuations in the system. The fluctuations are of course even more 
pronounced in the smaller N system. 

The random roughening of the interface makes the identification of the point of onset of the instability in the 
stochastic system somewhat subtle, given the small size of the pattern near onset. In order to compare the stochastic 
system to the MF prediction, we need a way to distinguish this random roughening of the interface from the gen- 
uine pattern forming mechanism discussed in the previous section. We devised two tests whose results we believe 
demonstrate the existence of an instability in the stochastic system. 

Both of these tests exploit the predicted transition between stable and unstable states that occurs when the channel 
width exceeds a critical value b*, as stated in (6). First, we measured the ensemble averaged velocity of the mean 
interface z = ^^ Nb ^^ ^ function of b (Fig. ^J. The increasing trend along the envelope of the different curves 
can be understood as a result of wider interfaces presenting an effectively larger number of particles N^ff. In fact, 
since in the steady-state, v ^ (IuNY^^, we see from the figure the remarkably simple result Ne/f ^ Nb. This simple 
dependence continues until b approaches b* (dotted vertical lines), where the velocity suddenly increases. This increase 
can be understood as a result of the system spending much of its time in a configuration in which one side of the 
interface significantly leads in front of the other. The lagging side then effectively stalls while the leading side is in a 
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FIG. 4: (Color online) Evidence of the transition to instability. The linear envelope of the different curves demonstrates the A*' 
renormalization Neff ~ Nb required by widening the channel. The dotted lines show the 6*resulting from the corresponding 
deterministic, MF cutoff simulation (fc=.25). Each data point represents an ensemble average of three long trials (100 time 
units each). The parameters are D = l, r-o = 6, a = 0.3. 



region of large reaction rate, and thus propagates quickly. The overall effect is an increase in the velocity averaged 
over the width of the channel. The fact that the change occurs so near the b* calculated earlier suggests that the 
cutoff approach is effectively capturing the stochastic dynamics. 

As a second test, we plotted the mean roughness of the interface H^ vs 6 (figure 3). M^ is defined in the standard 

way 
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{[{EjMhj)/N 
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where () denotes ensemble average and the bar denotes average over the transverse direction. For & < 6* we see power 
law scaling reminiscent of that discovered by Kardar, Parisi, and Zhang ^] for a growing interface. However, the data 
shows no sign of a universal exponent. It may be that the very weak stability of the interface near b* is responsible 
for a long crossover. This issue clearly requires more extensive study. For a fixed 6, W decreases with increasing TV, 
consistent with the hypothesis that interface roughness is noise driven in this regime. However, for b ^ b* this simple 
dependence is lost. The curves converge near b* , showing that particle number and its associated noise are no longer 
the relevant factor in determining interface roughness. Past this intersection, there is no apparent correlation between 
W and N. We interpret this as a crossover from noise driven interface roughness to gradient driven pattern formation 
occurring very near the b* predicted from the cutoff MF approach. 

Thus, the cutoff MF approach is quantitatively successful in predicting the transition of the width and velocity 
observed in ensemble averaged stochastic fronts when the channel is widened. In contrast, the naive MF approach 
predicts no such transition and an infinite steady-state velocity, in qualitative disagreement with the simulation results. 
The cutoff MF approach also predicts the velocity of the average interface for b < b*, provided we take N -^ Nb. 

Another aspect of the stochastic system which one would like to predict is the ensemble-averaged shape. We find 
that qualitatively this behaves as expected; namely, for small channel width the average shape is flat, and above 
the critical width, a nontrivial shape is apparent. The amplitude of the averaged pattern continues to increase with 
increasing width. However, we do not know how to quantitatively relate the average pattern to the results of the 
cutoff MF equations. One obvious impediment is the fact that the stochastic system switches parity at random, with 
the right and left sides alternating as the leading edge. Thus, a naive ensemble average produces a shape which is 
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FIG. 5: (Color online)Noise driven roughness scaling for b < b* , gradient driven scaling thereafter. Each data point represents 
an ensemble average of three long trials (100 time units each). Parameters are D = 1, ro = 6, a = 0.3. 



highest in the center, clearly at odds with the deterministic calculation. Another aspect of the problem that we would 
like to correlate with the deterministic calculation is the growth rate of the pattern near onset. This problem is also 
difficult because the width, measured in the usual way, consists of a contribution from noise driven roughening and 
one from pattern formation. Clearly, the usual MF approach can only make predictions about the contribution from 
patterning and thus the stochastic results and MF predictions are intrinsically difficult to compare. Even though 
the noise decreases with increasing TV, the dominance of the dynamics by the leading edge where fluctuations are 
unavoidably present makes this a nontrivial task, even at large N. These questions remain challenges for the future. 
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APPENDIX A: DETAILS CONCERNING NUMERICAL INTEGRATION 



In order to determine the spectrum Qo{v) we numerically integrated Eq. 5, and this required numerically integrating 
Eq. 3. Integration of Eq. 3 was initialized in the bulk state, from the left, where we defined z = 0. We took set v 
arbitrarily, took 0o(O) = .99, and calculated (J>q{0) from the solution to the version of Eq. 3 linearized around 00 ~ 1 
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■ Ol-u 
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T^init, which we set to one, differs from the previously defined rg in that it fixed the rate at a fixed location in the 
bulk state rather than where 0o = 1/2. Integration terminated in the neighborhood of the cutoff, half way between 
timesteps where (J)'q/(J)o crosses —v/D. N was then read off from the relation N = — j, , "",^ — r and the value of the 
cutoff was recorded for subsequent numerics. This value is measured relative to the bulk position where (J)q — .99, not 
relative to z, where (j>o — 1/2. All of this was done using ode45 from MATLAB, with a maximum stepsize of .001. 



The solution for 0o appears as a coefficient in Eq. (O , and was incorporated into the ODE integration scheme with 
a cubic sphne. Numerical integration of Eq. (jSJ was initialized at z = with some trial Qq, 77(0) = 1 and 



which follows from Eq. (0 if we plug in 0o ~ 1- Integration terminated at the Zcut, where we checked if Eq. 1)13(1 was 
satisfied with the trial fl. This procedure was iterated with a root solver while varying il until Eq. 1(13(1 was satisfied. 
Each integration of Eq. (jSJ was done over 1000 timesteps with a fourth order Runge-Kutta ODE solver with fixed 
step size, meant to facilitate incorporation of the spline. 

This yielded the exact numeric solution presented in Fig. 1. Our "exact analytic approximation" presented in Fig. 
1 is just the solution to Eq. ll3l obtained with a root solver. The required values for N and Zcut were obtained from the 
previous integration of Eq. , and we took tq = 1 . Since we dropped the term involving (pQ in Eq. |H1 this equation 
is insensitive to the precise definition of z, e.g. it could be defined naturally as the coordinate where (f>o — 1/2 or it 
could be defined out of numerical convenience as the coordinate where 00 = .99. This insensitivity explains why the 
results agree so well despite the fact that rg and rinu are not defined in the same way. 
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